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Abstract 

We consider three-dimensional convection of an incompressible fluid saturated in 
a parallelepiped with a porous medium. A mimetic finite-difference scheme for the 
Darcy convection problem in the primitive variables is developed. It consists of 
staggered nonuniform grids with five types of nodes, differencing and averaging op- 
erators on a two-nodes stencil. The nonlinear terms are approximated using special 
schemes. Two problems with different boundary conditions are considered to study 
scenarios of instability of the state of rest. Branching off of a continuous family of 
steady states was detected for the problem with zero heat fluxes on two opposite 
lateral planes. 
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Introduction 



Several algorithms ( finite element, finite difference, finite volume and spectral 
methods) are used for the simulation of physical problems involving partial 
differential equations (PDEs). Usually, the PDEs express fundamental physi- 
cal laws like conservation of mass, momentum and total energy. While solving 
such problems, some information about the problem and its structure is lost 
during the discretization that replaces the PDE by a system of algebraic equa- 
tions. In recent years so-called mimetic discretizations were developed to yield 
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stable and accurate solutions by preserving the analytical properties of the 
underlying PDEs [Tf2] . Several mimetic or conservative finite difference meth- 
ods were derived and their conservation properties were discussed for PDEs 
arising in fluid dynamics [3i4 ; 5j. Experience shows that discrete conservation 
of mass, momentum and kinetic energy produce better results as compared 
with nonconservative methods. 

The main goal of this paper is to develop a mimetic scheme for three-dimensional 
equations of the convection in a porous medium. Natural convection of incom- 
pressible fluid in a porous medium differs from the convection of a single fluid 
[6]. Usually, after the state of rest loses its stability, a finite number of regimes 
(convective patterns) may appear. An exciting example with an infinite num- 
ber of steady states was found for the planar problem of incompressible fluid 
convection in a porous medium based on the Darcy law [7] . This case of the ap- 
pearance of a continuous family of equilibria was explained by the cosymmetry 
theory [8ll9] . 

To compute a continuous family of steady states we need the numerical scheme 
to be mimetic of the underlying system. The approximation of the planar 
Darcy convection equations was done for the first time using the Galerkin 
method [10] . Then the computations of the families of steady states were per- 
formed using the finite-difference scheme [TTlfT2"] and a combination of spectral 
and finite-difference approaches [13] . In [13] a staggered grid discretization for 
the planar problem of Darcy convection was developed in primitive variables. 
Special approximation of the nonlinear terms of the underlying system was 
the crucial step in all these works. It was found that the loss of gyroscopic 
and cosymmetric properties in the finite-dimensional approximation destroys 
the family of steady states and leads to a finite number of isolated stationary 
regimes [TT][r5] . 

In this work we consider the three-dimensional problem of natural convection 
in a porous medium and develop a finite-difference scheme for a system in 
primitive variables (velocities, pressure and temperature). The discretization 
is based on nonuniform staggered grids [H] with five types of nodes by using 
the differencing and averaging operators on two-nodes stencil. The present 
finite difference scheme is constructed in the spirit of the fully conservative 
second-order finite difference scheme for incompressible flow on nonuniform 
grids [5J. An algorithm for the computation of the family of steady states is 
described, and is used in numerical computations. 

This paper is organized as follows. The equations for the three-dimensional 
Darcy convection problem in primitive variables are formulated in Section 
1. In Section 2, the grids, discrete operators and discretization in space are 
described. The computation of the continuous family of steady states with 
varying stability spectra (cosymmetric family) is described in Section 3. Nu- 
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merical results on the branching off of the family of steady states and isolated 
regimes from the state of rest are presented in Section 4. 



1 Darcy convection 

1.1 Darcy equations in primitive variables 

We consider an enclosure filled by a porous medium saturated by an incom- 
pressible fluid which is heated from below. We assume that the fluid is in- 
compressible according to the Boussinesq approximation The system of 
equations consists of the momentum equation based on the Darcy law 



-v = —Vp + #7 — v , 

e 

and the continuity and energy equations 
V ■ v = 0, 

= A6 - Xv • 7 - (v ■ V)0. 



(1) 

(2) 
(3) 



Here v = (v 1 , v 2 , v 3 ) T is the velocity vector, 7 = (0, 0, — 1) T is the direction of 
gravity, p(x, y, z, t) is the pressure, 8(x, y, z, t) is the deviation of the temper- 
ature from a linear (in z) profile, e is the porosity of the medium, x, y, z are 
the space variables and the dot accent denotes differentiation w.r.t. time t. 

The Rayleigh number is defined as A = agTKl/x^, where a is the thermal 
expansion coefficient, g is the gravity acceleration, \x is the kinematic viscosity, 
X is the thermal diffusivity of the fluid, K is the permeability coefficient, T 
is the characteristic temperature difference and I is the length parameter. We 
suppose that the temperature at the boundary is given by a linear function 
on the vertical coordinate z. 

The parallelepiped V = [0,L X ] x [0, L y ] x [0, L z ], with length L x , depth L y 
and height L z is filled with the fluid. The normal component of the velocity 
is equal to zero at the boundary 

V-n = 0, (x, y, z) G dV. (4) 

We consider two problems with different boundary conditions for the temper- 
ature. The first problem is characterized by a temperature deviation 9 equal 



to zero at the boundary &D (Dirichlet boundary condition): 
9 = 0, (x, y, z) G &D. 



(5) 



The second problem has mixed boundary conditions: the heat flux equals zero 
on two lateral faces d\D = {y = 0} U {y = L y } and the temperature deviation 
9 is equal to zero on the remaining faces d 2 D = D \ d\D 

9 y = 0, {x,y,z)ed 1 D, 9 = 0, (x,y,z)ed 2 D. (6) 

The initial condition is given as follows 

9(x,y,z,0) = 9 (x,y, z), v(x, y, z, 0) = v (x, y, z). (7) 

It is simple to check that the equations ([I])-© are invariant with respect to 
the discrete symmetries 

Rx ■ W 3 ,p,6>} i-> {L x - x,y,z,-v l ,v 2 ,v 3 ,p,9}, (8) 

R y : {x,?/,2:,t;W 3 ,p,6)} i-> {x,L y - y, z, v \ -v 2 , v 3 ,p, 9}, (9) 
Rz ■ {x,y, z,^ 1 ,^ 2 ,^ 3 ,^,^} ^ {x,y,L z - z,v\v 2 ,-v 3 ,p,-9}. (10) 

This implies that with appropriate transformations of velocity, pressure and 
temperature deviation, a given set of solutions produces a new set. 

1.2 Darcy equations for the temperature and stream function 

When the initial velocity v q and the initial temperature distribution 9 do not 
depend on y the system ([I])-©, © has a two-dimensional solution 

v 1 = v 1 {x, z, t), v 2 = 0, v 3 = v 3 (x, z,t), p = p(x, z,t), 9 = 9(x, z, t). (11) 

In this case we can write our system as a system containing temperature and 
stream function. We follow the usual assumption in porous media flow and 
neglect inertia in the momentum equation. The continuity equation ([2]) is 
fulfilled automatically when the stream function ip is given by 

v 1 = -i/> 9 , v 3 = if> w . (12) 

Then the underlying system can be transformed to a new form. After appli- 
cation of the cur/-operator to ([T|) we deduce 

= A 2 if>-9 x = G, A 2 Q = Q ts + Q gg , (13) 
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and using (fT2j) we obtain from ([3]) 



e = A 2 e + \i, x + j(e, v) = f, j(0, ^) = e x i> g - e z ^ x . (14) 



The boundary conditions for the system (IT5I) . ( THl) follow from (J6j): 

^ = ^ = 0, (a;, 2) G 9P, where X> = [0, L x ) x [0, Lj. (15) 



and the initial condition is formulated only for the temperature 

e(x,z,o) = e (x,z), (16) 



where 6q denotes the initial temperature distribution. For a given 6q, the 
stream function ip can be obtained from (fT3l) and ([6]) as the solution of the 
Dirichlet problem via Green's operator ip = G9 X . 

Equations (fT5j) - (}T5|) require that the equilibrium 6 = ip = (state of rest), 
be stable if A < An, where A nm = (2-nnja) 2 + (2irm/b) 2 (m,n e Z) are the 
eigenvalues for the corresponding spectral problem. It was shown in [9] that 
the first critical value An has multiplicity two for any domain D. As a result, 
a continuous family of steady states appears [7115] . 

The system f lT3|) - f|T5l) possesses the cosymmetry property [8]: a vector-function 
(9, —ip) being orthogonal to the right-hand side of ( |T3|) and < HM in L 2 . Then, 
we obtain the cosymmetry condition in the following form 

J (F^ - G6)dxdz = J [A6ip - Aip6 + Xip x ip + 6 X 6 + J(6, xp)ip] dxdz = 0.(17) 
v v 



This can be checked directly using integration by parts and Green's formulae. 

For the integration of equations f|T3|) - f|T5|) it is essential to provide a discrete 
version of the cosymmetry condition. In [11] a regular uniform mesh was used 
and both stream function and temperature were defined at the same nodes. 
The Jacobian approximation was based on the Arakawa scheme and a 
number of one-parameter families of steady states were computed. It was also 
found that a violation of the cosymmetry property led to a degeneration of 
the family. The application of staggered nonuniform grids for the problem 
( fT3l) - (|T6l) was considered in [T3] . 
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2 Spatial discretization 

We have discretized the equations (JTJ) — (ED using five different types of nodes: 
one for the pressure, another for the temperature and three nodes for the 
components of the velocity vector, see Fig. HJ 



z 



X 




Fig. 1. Grid and nodes 

For the first problem with Dirichlet boundary conditions (jl]) and §5$) we intro- 
duce a nonuniform grid for the temperature 9 

= x < xi < . . . < x Nx+1 = L x , 
= y <y 1 < ... < y Ny+ i = L y , 
= z < z 1 < ... < z Nz+1 = L z . 

The grid for the second problem with mixed boundary conditions (jl]) and 
differs only in y direction 

2/o + 2/i = 0, yi< ... <y Ny , y Ny + y Ny+1 = 2L y . 

Thus the temperature 9 is defined at the nodes 

u = {(xi, yj, z k ), i = 0, . . . , N x + 1, j = 0, . . . , N y + 1, k = 0, . . . , N z + 1}. 

We introduce then the staggered grids along all coordinates 
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Vj+1/2 



%i+l/2 




Zk+1/2 = -j{%k + Zk+l), k = 0,...,N x . 

The velocities v l , v 2 and v 3 are defined on the grids which are staggered 
respectively along the corresponding coordinates 

wi = {(xi,y j+1/ 2,z k+1/2 ), i = 0, . . . , N x + 1, j = 0, . . . , N y , k = 0, . . .,N Z }, 
^2 = {(x i+1/ 2,yj,z k+1/2 ), i = 0, . . ., N x , j = 0, . . . , N y + 1, k = 0, . . . , N z }, 
^3 = {(^+i/2 ; %-+i/2,^), i = 0, . . . ,N X , j = 0, . . . , N y , k = 0, . . . , N z + 1}. 

Finally, the pressure p is defined at the nodes 

= {{ x i+i/2, V Zk+1/2), i = 0, . . . , N x , j = 0, . . . , Ny, k = 0, . . . , N z }. 

In the case of Dirichlet boundary conditions we do not need any fictitious 
nodes for the discretization. In the case of mixed boundary conditions the 
grids are introduced in such a way that on d%D the boundary conditions for the 
temperature and the normal component of velocity are fulfilled automatically 
We define fictitious nodes for the temperature and velocity v 2 to approximate 
the boundary conditions on d\D. 

2.1 Discrete finite-difference operators 

To approximate (JTJ) — (EJ) we define a set of discrete analogs of first order differ- 
ential operators on a two-point stencil 



dl9i+l/2,j,k 



@i+l,j,k @i,j,k 




x)i+l/2,j,k, 



d2@i,j+l/2,k 



Vj+i - Vj 



(Qy)i,j+l/2,k, 



(18) 



d^i,j,k+l/2 




{&z)i,j,k+l/2, 



Zk+l — z k 



and weighted averaging operators on the coordinates 



$l9i+l/2,j,k 



[%i+l — Xi+l/2)9i+l,j,k + {%i+l/2 ~ %i)9i,j,k 
%i+l 



(Q)i+l/2,j,k, 
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(s/j+i ~ yj+i/2)0i,j+i,k + {yj+i/2 - yj)Qi,j,k _ /n\ /iq\ 

Vj+i - Vi 

x D _ ( z k+l — z k+l/2)9i,j.k+l + { z k+l/2 — z k)9i,j,k , p s 
<J3^i,j,k+l/2 — ~ [V)i,j,k+l/2- 

The formulas (fT8l) - (|T9l) are valid both for integer and half- integer values of i, 
j and k. Then the discrete analog of the Laplacian on the seven-nodes stencil 
can be written as 

A h = didi + d 2 d 2 + d 3 d 3 « A, (20) 



and the averaging operator on three-dimensional cell is given as 

So = S 1 5 2 5 3 . (21) 



The nonlinear term approximation is constructed using a linear combination 
of two terms 



{vV0) itj , k nJ{6,v) iijtk = (22) 
' \dM95 2 5 3 v l ) + d^OS^v 2 ) + d 3 S 3 (9S 1 5 2 v 3 ) 



+ 3 



2 r 

d x b 2 b 3 {b$hv x ) + d 2 5i5 3 (5o^ 2 v 2 ) + d^^oOdzv 3 ' 



This provides second order accuracy for the uniform grid and an asymptoti- 
cally second order accuracy for the quasi-uniform mesh. It allows conservation 
of energy and constitute a mimetic discretization of the underlying problem. 



2.2 Semi- discretization 



To reach some steady state it is useful to apply an approach based on artificial 
compressibility [To^TIT] . Thus instead of equation (J2J) we consider the following 
equation with coefficient rj 

p+iv-w = 0. (23) 
V 

Using the operators (fT8l)- (l22l we discretize system (CEJ), (j3J) and (1231 in the 
following form 



6-A h 6- XSiS 2 v 3 + J{9, v) = 0, (24) 

J i,j,k 
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v 1 — e(—dip — v 1 ) 
v 2 — e(—d,2P — v 2 ) 



i,j+l/2,fc+l/2 



i+l/2J,fc+l/2 



v 3 - e(-d 3 p - V 3 + 5^20) 

p + -(d l v l + d 2 v 2 + d 3 v 3 ) 
V 



i+l/2J+l/2,fc 



J i+l/2j'+l/2,fc+l/2 



0. 



(25) 
(26) 
(27) 

(28) 



For the problem with Dirichlet boundary conditions ([I])-© the grids are in- 
troduced in such a way that the nodes for the transversal velocity are located 
on the wall. It allows us to fulfill the boundary conditions on a rigid wall (jSj) 
in a very simple way for each pair of planes: 
for x = (i = 0) and x = L x (i — N x + 1) we have 



v 



l 

i,j+l/2,fc+l/2 



0, j = 0,...,N v , k = 0,...,N z , 
di,j,k = 0, j = 0, . . . , N y + 1, k = 0,...,N z + l, 

for y — (j = 0) and y = L y (j — N y + 1 )we have 



(29) 



i+l/2j',fc+l/2 



0, i = 0,...,N x , k = 0,...,N„ 



0i,j,k = 0, z = 0, . . . , N x + 1, k = 0,...,N z + l, 
and for z = (k = 0) and z = L z (k — N z + 1) we have 



(30) 



i+l/2j'+l/2,fc ~~ u ' 



0, i = 



A^, j = 0, iV„ 



0, i = 0, 



j • • • > 

A^ + l, j = 0,...,JV w + l. 



(31) 



The problem with mixed boundary conditions (JXJ) — (J4j) and ([6]) is discretized 
using fictitious nodes to satisfy the boundary conditions on the planes y = 
and y = L y . Thus instead of (13"01) we have 

^+1/2,0^+1/2 = — u £t-i/2,i,fc+i/2> i — 0,..., N x , k — 0, . . . , iV^, (32) 

U i+l/2,JV v +l,Jfe+l/2 = ~~ u i+l/2,iV H ,fe+l/25 i = 0, . . . , N x , k — 0, . . . , iVg, 

0i,O,fc = ^i,l,fc> 9i,N y +l,k — 9i,N y ,k, i = 0, . . . , N x + 1, = 0, . . . , N z + 1. 

Discrete equations for the planar problem 



The system of ordinary differential equations ([24]) -([28]), (|29l), (|3T]) and (|32l) 
has a solution such that ^+1/2^,^+1/2 = an d the other variables are invariant 
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w.r.t. index j. Then, we can exclude equation (126!) and consider the two- 
dimensional problem. Moreover, we can eliminate pressure and velocities and 
reformulate the system with the discrete stream function and temperature. 
Let's define the stream function if) at the nodes u> using difference operators 

tfHD 

= -vJfc+i/2. {diip)i+i/2,k = vf +1/2tk . (33) 



The discrete analog of the continuity equation (j2J) is automatically fulfilled 
by ( 1331) . After the substitution of ( 1331) in ( 125]) and ( 1271) and the combination 
of the resulting formulas we can deduce the analog of ( fTBl) (inertia terms are 
omitted) 

[A a>h V - A^]^ = 0. (34) 
Similarly we find from 



i k 



1 2 
A 2 , h # + XD x iP + -J D + -J d 



i,k 



(35) 



Here D x = d x 8i, D z = d 3 8 3 are the first order differencing operators on three- 
nodes stencils. The Laplacian and Jacobian are approximated using 



A 2 , h = dtdi + d 3 d 3 , 

J D = D x (9D y iP) - D y (9D x iP) , (36) 
J d = d x (d Bd z if^ - d z (doBd x if)J , 

where d = 5i5 3 , d x = d x S 3 , d z = d z 5 x . 

The resulting scheme (I34l) - fl36l) ensures the fulfillment of a discrete analog of 
the cosymmetry property (ITT)) as well as the nullification of the gyroscopic 
terms |12j . Equations in [TTJ follow from (l3Tj) - (l3"6)) for the case of the uniform 
grid Xi = ih, z k = kg, h = L x j (N x + 1), g — L z j (N z + 1). Then the Jacobian 
approximation gives the famous Arakawa formula [15J on uniform grids h — g. 



3 Computation of the family of steady states 



We rewrite the resulting system of equations in vector form. Let us introduce 
vectors which contain only unknowns at internal nodes 
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e = ( 


,^1,1,1' • • • ' 


9n x ,i,i, 6*1,2,1, • • • , 9n x ,n v ,n z ) 


• 


v l = 


(^1, 1/2, 1/2: 


• • • > v N x ,l/2,l/2i ^1,3/2,1/2; ■ ■ ■ 


5^,^+1/2,7^+1/2)5 


v 2 = 


(^1/2, 1,1/2: 


2 2 
■ ■ ■ ' ^^+1/2,1,1/2 5 V l/2,2, 1/25 • 


• • 5 u N x +l/2,N y ,N z +l/2 


v 3 = 


(^1/2, 1/2,1: 


3 3 
' • • • 5 v N x +l/2,l/2,l' ^1/2,3/2,15 • 


• • 5 u N x +l/2,N y +l/2,N z 



P — (Pl/2,1/2,1/2, • • • 5^+1/2, 1/2,1/2 5 Pl/2,3/2,1/2 5 ■ • ■ , PN x +l/2,N y +l/2,N z +l/2) , 

and obtain the system which corresponds (!24{) - (l28l) 

= A 1 Q + XC 1 V 3 -F{Q,V), 
V 1 = -B 4 P-C 2 V 1 , 
Vi = -B 5 P-C 3 V 2 , 

V 3 = -B 6 P - C A V 3 + C 5 0, (37) 
P = —B1V 1 - B 2 V 2 - B 3 V 3 . 

Here the matrices Bf., k = 1, . . . , 6, are constructed by the application of first 
order difference operators, and the matrices Ck, k = 1, . . . , 5, by the averaging 
operators. The matrix A\ presents the discrete form of the Laplacian. The 
nonlinear term is given by -F(0, V). Equations (|37|) form a system of 

5N x N y N z + 3(N x N y + N X N Z + N y + N z ) + 2{N X + N y + N z ) + 1 

unknowns. 

From (1371) at J = we can derive the perturbation equations (a is a decrement 
of linear growth) to analyze the stability of the state of rest 



a0 = Ai6 + XdV 3 , (38) 

aV 1 = -B 4 P - C 2 V\ (39) 

aV 2 = -B 5 P - C 3 V 2 , (40) 

aV 3 = -B 6 P - C 4 V 3 + C 5 0, (41) 

aP = -B x V l - B 2 V 2 - B 3 V 3 . (42) 



For the decrement a = we obtain the system from which we can determine 
the threshold value of the Rayleigh number corresponding to the monotonic 
loss of stability. We can express P, V 1 , V 2 V 3 via 6 from fl3"gl)-fT4"2D and 
obtain a system of N x N y N z equations for the unkown vector ©. Substituting 
(I3"9l, P0|) and (jill into (@2]) we deduce 

A 1 e = XC 1 (C 5 -B 6 S)Q. (43) 
Here we find the vector P = SO from the system of linear algebraic equations 
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with rank deficiency 1 

{B\C^ B± + B 2 C 3 X B 5 + -B3C 4 1 B 6 )P = 5 3 C 4 r 1 C 5 0. 

Since for an incompressible flow the pressure may differ by a constant we can 
exclude one component of P and respectively one equation. 

To find an isolated convective pattern we apply the direct approach and inte- 
grate the system of ordinary differential equations (EJTI) by the classical fourth 
order Runge-Kutta method up to convergence. 

To compute a family of steady states we apply the technique based on the 
cosymmetric version of the implicit function theorem [18] and the algorithm 
developed in [T0llll|ll9j . The zero equilibrium V 1 = V 2 = V 3 = P = = 
is globally stable for A < Ai where Ai is the minimal eigenvalue for spectral 
problem (I43p . When A is slightly larger than Ai, then all points of the fam- 
ily are stable [5]. Starting from the vicinity of unstable zero equilibrium we 
integrate the ordinary differential equations fl37|) up to a point 0o close to 
some stable equilibrium on the family. Then we correct the point ©o by the 
Newton method. To predict the next point on the family we determine the 
kernel of the linearization matrix (Jacobi matrix) at the point Go and then 
use the Adams-Bashford method. This procedure is repeated to obtain the 
entire family of steady states. It is important to note that the given procedure 
allows us to compute the stable regimes as well as unstable ones. 



4 Numerical results 

V. Yudovich [5] proved that the appearance of a family of steady states in 
the planar Darcy convection is caused by the nontrivial cosymmetry of the 
problem. Loss of stability of the state of rest is characterized by the repeated 
eigenvalues for the corresponding spectral problem. For the problem under 
consideration we find the critical Rayleigh numbers from the system (l4"3"j) . 
There exist two scenarios of instability of the state of rest in the parallelepiped: 
branching off of isolated regimes and the appearance of a family of steady 
states [20]. In our computer experiment the emergence of the family was only 
observed for the problem with mixed boundary conditions. It was found that 
a family of stable equilibria has appeared in the case of rather small value of 
L y /L x (relative depth) which depends also on the ratio L x /L z . 

It was shown in [14j that a uniform grid is the best choice for the computation 
of the critical Rayleigh numbers while a nonuniform grid is useful for the 
computation of a specific regime with a desirable accuracy. Thus, we use the 
uniform grids and set up the amount of internal nodes for the temperature as 
N x x N y x N z . 
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4-1 Critical Rayleigh numbers 



We present in Table [T] the first seven critical Rayleigh numbers A for the 
problem with Dirichlet boundary conditions for several values of the depth 
L y and fixed length L x = 2 and height L z = 1. In the case L y < 1 [L y > 1) 
the computation of critical Rayleigh numbers was carried out on the mesh of 
14 x 6 x 6 (14 x 10 x 6) internal nodes for the temperature. 

Table 1 

Dependence of critical Rayleigh numbers on the depth L y for the problem with 
Dirichlet boundary conditions; L x = 2, L z = 1, mesh 14 x 6 x 6 (14 x 10 x 6) 



Ly 


0.4 


0.6 


0.8 


1.2 


1.6 


Ai 


157.5 


99.2 


77.9 


62.0 (59.1) 


(53.5) 


A 2 


159.0 


100.1 


78.0 


62.3 (59.7) 


(54.9) 


A 3 


207.0 


138.4 


108.5 


83.4 (73.3) 


(60.6) 


A 4 


209.4 


139.6 


116.0 


85.1 (78.2) 


(66.0) 


A 5 


298.3 


178.7 


122.9 


90.0 (79.6) 


(69.5) 


A 6 


303.8 


191.1 


132.4 


100.1 (95.2) 


(81.4) 


A 7 


310.5 


192.0 


152.1 


107.6 (95.3) 


(86.9) 



One can see that some critical values in Table [T] are close to each other. When 
we take L x = L y we find that some critical values coincide. For example, for 
D = [0, 2] x [0, 2] x [0, 1] and the mesh 12 x 12 x 6 we find that Ai = 51.3 
and A2 = A3 = 53.8. But this is a consequence of the discrete symmetries on 
the Xi and x 2 coordinates and doesn't lead to the appearance of a family of 
steady states. 

We summarize in Table [2] the first seven critical Rayleigh numbers A for the 
mixed boundary conditions problem for several values of the depth L y and 
fixed length L x = 2 and height L z = 1. One can see that for L y = 0.4 
and L y = 0.6 two minimal eigenvalues of the problem (l4~3!) are repeated. 
This corresponds to the birth of the continuous family of steady states. The 
parallelepiped with L y = 0.8 gives a single minimal eigenvalue and that results 
in the appearance of two isolated stationary regimes. 
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Table 2 

Dependence of critical Rayleigh numbers on the depth L y for the problem with 
mixed boundary conditions; L x = 2, L z = 1, mesh 14 x 6 x 6 (14 x 10 x 6) 



y 


0.4 


0.6 


0.8 


1.2 


1.6 


Ai 


52.5 


52.5 


46.7 


44.5 (42.5) 


(45.7) 


A 2 


52.5 


52.5 


52.5 


51.6 (49.4) 


(46.1) 


As 


90.8 


56.6 


52.5 


52.5 (52.5) 


(49.6) 


A 4 


93.1 


66.8 


56.2 


52.5 (52.5) 


(52.5) 


A 5 


93.1 


84.7 


73.2 


67.8 (58.1) 


(52.5) 


A 6 


102.1 


93.1 


93.1 


68.3 (64.8) 


(57.8) 


A 7 


122.2 


93.1 


93.1 


80.9 (68.6) 


(66.3) 



It should be noted that on the fixed mesh the threshold corresponding to the 
branching off of the family of steady states (repeated critical values) doesn't 
depend on the depth L y . Comparison with results on the mesh 24 x 6 x 12 
shows that even a rough mesh allows to find a threshold (minimum of critical 
Rayleigh numbers) with accuracy about 10%. For instance, using the mesh 
24 x 6 x 12 we obtain A x = 50.5 for L y = 0.4 and L y = 0.6 and X 1 = 46.8 for 
L y = 0.8. 

On the other hand convection in a three-dimensional box with insulating im- 
permeable lateral boundaries (the lateral walls are taken as thermally insulat- 
ing) was a subject of numerous works [6]. We may refer here to [21] where it 
was shown that different two-dimensional and three-dimensional states appear 
depending on the initial conditions. 



4-2 Computation of steady states 



Now we study the mixed boundary conditions problem and consider the par- 
allelepipeds with small depth when the family of stationary solutions branches 
off. The steady states belonging to the family are essentially two-dimensional 
and don't depend on the coordinate y. 

Fig. [2] demonstrates this by the presentation of several convective patterns 
from the family: none has transversal motion to the plane y = const. The 
relative location of each steady state is given in Fig. [3] where each letter cor- 
responds to a flow pattern in Fig. [2] In Fig. [3] we use the Nusselt values for 
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the planar problem [TO] 

Nu v = J 9 x (^,0,z)dz, Nu h = J 9 z (x,0,0)dx. (44) 



Here Nu v corresponds to the cumulative heat flux from left to right defined 
for the centered vertical section of the rectangular domain. The value Nuh is 
a combined heat flux through the bottom of the enclosure, y = 0. 








Fig. 2. Members of the family of steady states; A = 60, D = [0, 2] x [0, 0.5] x [0, 1] 

The regimes in Fig. 2 are the members of the family of steady states, because 
stability spectra for each state have exactly one value being zero with reliable 
accuracy (1CT 8 ). We plot the distribution of eigenvalues of the Jacobi matrix 
(matrix of linearization) computed for the convective pattern with two sym- 
metrical rolls (regime a in Fig. [3]). We take different values of the depth L y 
to demonstrate the three-dimensional instability on the family, see Fig. H]). 
It is clearly seen that exactly one point a lies on the imaginary axis. The 
corresponding eigenvector defines the neutral direction along the family. Such 
a family is called cosymmetric, the stability of its members is governed by 
nonzero eigenvalues or defined on the submanifold being transversal to the 
family. 

It should be noted that the given Rayleigh number is rather far from the 
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Fig. 3. Family of steady states; A = 60, D = [0, 2] x [0, 0.5] x [0, 1] 
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Fig. 4. Stability spectrum for the two-rolls steady state from the family; A = 100, 
L x = 2, L z = 1 

critical value when some members of the family become unstable. In the case 
of the parallelepiped D = [0, 2] x [0, 0.5] x [0, 1] the family of steady states 
appears at A ~ 51. The transformation of stability spectra with increasing A 
is presented in Fig. [5J One can see that the symmetrical convective pattern 
lost its stability before A = 200. Some states on the family become unstable 
at A ~ 190. This value is less than A ~ 400 which corresponds to the critical 
value of instability for the planar problem. 



We present in Fig. [6] the distribution of spectra for different values of the 
depth L y . One can see that the steady state with two rolls is stable when 
L y = 0.3 and unstable for L y > 0.4 (two spectra values in right part of the 
complex plane). When the Rayleigh number becomes greater the instability 
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Fig. 5. Stability spectrum for symmetric two-rolls pattern from the family for dif- 
ferent A; D = [0, 2] x [0, 0.5] x [0, 1] 
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Fig. 6. Stability spectrum for the two-rolls steady state from the family; A = 200, 
L x = 2, L z = 1 



When the depth L y = 0.9 the only isolated stable regimes are branched off of 
the state of rest. Because of the discrete symmetries of the problem fl8|)- f|T0|) we 
obtained two regimes. The flow pattern for one of these regimes is presented in 
Fig. [3 The distribution of eigenvalues for this stable steady state is displayed 
in Fig. [HJ One can see that the distribution of the stability spectra has no 
point close to the imaginary axis. This regime is essentially three-dimensional 
and isolated. Branching off the isolated convective regimes is characteristic for 
the parallelepipeds with non-small depth. 
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Fig. 7. Isolated steady state; A = 52, D = [0, 2] x [0, 0.9] x [0, 1] 
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Re a 

Fig. 8. Stability spectrum for the isolated steady state from the family; A = 52, 
D = [0, 2] x [0, 0.9] x [0, 1] 

5 Conclusion 

Convection in a porous parallelepiped has demonstrated two different scenarios 
of instability of the state of rest. It was found by computer experiments that 
for zero heat fluxes on two lateral sides (mixed boundary conditions) the 
appearance of a cosymmetric continuous family of equilibria becomes possible. 
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To compute a continuous family of equilibria one needs to solve repeatedly a 
nonlinear system of algebraic equations that is degenerated in the vicinity of 
the family. This is why discretization is so important for the Darcy convection. 
We have developed the approach based on primitive variables equations and a 
finite-difference discretization with staggered nonuniform grids. This scheme 
mimics the nontrivial characteristics of the underlying problem that admits 
existence of a continuous family of steady states. 
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